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I. Introduction 

Simulation of multibody dynamics systems - such as robotic manipulators, automo- 
biles maneuvering and satellites deployment - remains a challenge to the dynamist due to 
its increasing roles in design improvements, control and safe operation. Because of sub- 
stantial progress made during the past three decades in formulation 1 - 19 , constraint treat- 
ment absolution techniques 21 " 36 and the availability of multibody dynamics simulation 
packages , it has now become almost a routine practice to perform realistic modeling 
and assessment of some practical problems such as mechanical linkages and manipulations 
of robotic arms if multibody components consist mostly of rigid bodies, discrete springs and 
dampers (see, e.g., Haug 15 ). However, substantial advances in modeling, formulation and 
computational methods are necessary in order to develop a real-time simulation capabil- 
ity for ground vehicle maneuvering dynamics, robotic manipulations and space structures 
deployment / assembly. 

Specifically, improved modeling of flexibility for localized motions and geometric non- 
linearities, material nonlinearities and contact /friction phenomena, robust and accurate 
treatment of the system constraint conditions and efficient use of emerging computer hard- 
ware/software technology continue to offer intense research opportunities. Thus, the de- 
velopment of a real-time multibody dynamics simulation capability requires a concerted 
integration of various modeling, formulation and computational aspects. These include: 
selection of a data structure for describing the system topology, computerized generation 
of the governing equations of motion, implementation of suitable solution algorithms, in- 
corporation of constraint conditions and easy interpretation of the simulation results. Of 
these, this chapter is concerned with three computational aspects of multibody dynamics 
simulation: direct time integration of the governing equations of motion, stabilization of 
constraint solution process and their computer implementation aspects. 

From the computational viewpoint, multibody dynamics (MBD) problems are distinct 
from the structural dynamics problems in that the solution of MBD problems must also 
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satisfy, at each time integration step, the attendant kinematic and equilibrium constraints. 
This has motivated many dynamists to develop various techniques, in addition to direct in- 
tegration algorithms, for accurately and efficiently handling the system constraints. Hence, 
reliability and cost of existing MBD simulation packages have been strongly affected by 
how efficiently and accurately the constraints are preserved during the numerical solution 

stage. 

In general, there have been two types of direct time integration algorithms for the 
transient response analysis of dynamical systems: explicit and implicit algorithms (see, 
e.g., Hughes and Belytschko 43 , Park 44 and Belytschko, Englemann and Liu 45 ). Currently, 
implicit algorithms appear to be favored by many MBD specialists when both the gen- 
eralized coordinates and the constraint forces axe treated as the unknowns. In this case, 
the corresponding formulations incorporate the system constraints by the Lagrange mul- 
tipliers method. It has been well known that the resulting Newton-like solution matrix is 
stiff. This has led to implicit time discretization of the constraint-augmented equations 
and simultaneous solution of both the generalized coordinates and the Lagrange multipli- 
ers. This approach has been extensively investigated by Gear 21 , Baumgarte 1 Orlandea, 
Chase and Calahan 23 , Petzold 27 , Nikravesh 31 , among others. Because these methods solve 
both the generalized coordinates and the constraint forces simultaneously, they will be 
called the simultaneous solution methods in this chapter. 

On the other hand, if the constraints axe eliminated so as to reduce the number of 
unknowns, it is possible for one to employ either implicit or explicit algorithm. For this 
situation, one may invoke either a geometric or algebraic procedure to streamline the re- 
sulting equations of motion if the system topology is an open tree. In essence, geometric 
procedures have utilized an open-tree topology such as the use of the incidence matrix by 
Wittenburg 10 and the body array matrix by Huston 19 . Some of the proposed algebraic 
procedures include the singular decomposition by Walton et al , the use of the general- 
ized speed of Kane and Levinson 20 , the coordinate partitioning technique by Wehage and 
Haug 28 , the selection of independent coordinates through the natural-coordinate formu- 
lation of Garcia de Jalon et al 33 and the so-called order-N procedures of Armstrong 11 , 
Holler bach 12 , Schwertassek and Roberson 17 , Orin, et al 25 , among others. 

As the complexity of MBD systems increases, the simultaneous solution methods 
have become less attractive. This is due to matrix ill-conditioning especially for the so- 
called index two and higher index problems (see, e.g., Ref. 27 and Brenan, Campbell 
and Petzold 46 for the definition of index for constraint characterization), divergence of the 
solution away from the constraint conditions, and ultimately, due to a large size of the 
equations that must be handled. As an alternative to the simultaneous solution methods, 
a series of computational methods that employ a divide- and- conquer strategy have been 
developed, which are termed as partitioned solution procedures presented in Park 47 , Felippa 
and Park 48 and Park and Felippa 49 . As an example, partitioned solution procedures allow 
one to analyze fluid-structure interaction problems with two separate single-field analysis 
packages, namely, the structural dynamics module and the fluid dynamics analyzer. At 
each time integration step, one may advance the solution of structural equations of motion 
by treating the fluid coupling term as an external force. Once the structural coordinates 
are advanced, the fluid state variables can be advanced by treating the structural coupling 
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terms as a source term. A naive partitioned procedure, however, can suffer from a loss of 
accuracy as well as computational stability. Thus, a combination of equation augmentation 
and stabilization should be devised to recover the accuracy loss and maintain unconditional 
stability. Such a solution procedure is in contrast to a practice of embedding both the 
structural and fluid dynamics attributes into a combined analysis program. 

The numerical solution procedure for MBD systems which we advocate in this chapter 
is ermed a staggered MBD solution porcedure that solves the generalized coordinates in a 
separate module from that for the constraint force. This requires a reformulation of the 
constraint conditions so that the constraint forces can also be integrated in time. A major 
advantage of such a partitioned solution procedure is that additional analysis capabilities 
such as active controller and design optimization modules can be easily interfaced without 

^follows 8 them mt ° a m ° nolithic P ro S ram - To this end, the rest of the chapter is organized 

After introducing the basic equations of motion for MBD system in the next sec- 
tion, Action III briefly reviews some constraint handling techniques and introduces the 

staggered stabilized technique 34 - 35 for the solution of the constraint forces as independent 
variables. 

. The numerical direct time integration of the equations of motion is described in Sec- 
tion IV. As accurate damping treatment is important for the dynamics of space structures 
we have employed the central difference method and the mid-point form of the trapezoidal 
rule since they engender no numerical damping. This is in contrast to the current prac- 
hce in dynamic simulations of ground vehicles by employing a set of backward difference 
formulas . First, the equations of motion is partitioned according to the translational and 
the rotational coordinates. This sets the stage for an efficient treatment of the rotational 
motions via the singularity-free Euler parameters. The resulting partitioned equations of 
motion are then integrated via a two-stage explicit stabilized algorithm for updating both 
the translational coordinates and angular velocities 34 . Once the angular velocities are ob- 

ained the angular orientations are updated via the mid-point implicit formula employing 
the Euler parameters. r J 6 

When the two algorithms, namely, the two-stage explicit algorithm for the generalized 
coordinates and the implicit staggered procedure for the constraint Lagrange multipliers, 
are brought together m a staggered manner, they constitute a staggered explicit-implicit 
procedure which are summarized in Section V. Section VI presents some example problems 

and discussions concerning several salient features of the staggered MBD solution procedure 
are ottered in Section VII. 


II. Governing Equations of Motion 

The Lagrangian equations of motion for mechanical systems that are free from any 
constraint can be written, for the generalized coordinate component u,, as 

c[dL dL 

— Qn i = l ... n. (1) 


dt diii dui 

where L is the system Lagrangian, t is the time, ( ' ) denotes time differentiation and O, 
.s the generalized applied force. It is well-known that, if there are m-constraint conditions 
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imposed on {uj, 


i = 1 . . . n}, the above equation must be modified as 




( 2 ) 


k-1 


where A is the Lagrange multiplier and B ki is the i-th gradient component of the fc-th 
constraint equation, viz, for configuration constraints 


d$ k 

*fc(u) = 0, = k=l...m 


( 3 ) 


and for motion constraints 


d$ k 


#*(u,u) = 0, B ki = -^r~, k — 1 . . .m. 


( 4 ) 


Therefore, regardless of the nature of constraints one may express the equations of 
motion with constraints in the following form: 


M B T 
B 0 


](;)-(?) 


(5) 


where M is a positive-definite matrix and c depends on the nature of constraints. For 
example, for configuration constraints we have 


d ( d$ d d$ 
c ~~du^du U ^ 2 dt^du U) dt 2 


and for motion constraints 


c = 


d$ 

' dt 


( 6 ) 


( 7 ) 


An implicit time integration formula to solve (5) may be written as 


{ 


u n = 8u n + K 
u" = 8u n + h" 


( 8 ) 


where 8 is a stepsize that is dependent on the choice of formula, and h? and h^are 
formula-dependent historical vectors that consist of past-step solution components - . 

As an example, the trapezoidal rule has the following 8 and historical vectors 


8 = h/2 

h" = ii "" 1 -Mu ”” 1 
h" = u- 1 + 6u n_1 


( 9 ) 


where h is the time-step increment. 
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Substitution of (8) into (5) yields 


'M S 2 B 1 "\ /u"\ / _ (Mh n + 6 2 Q n \ 

B 0 J V A ”/ V r A/ _ V Bh n + S 2 c ) 


( 10 ) 


In practice, in order to avoid pivoting and to maintain high accuracy, the solution of 
the above difference equations is carried out as follows. First, since M is nonsingular for 
properly formulated dynamical problems, one computes 


u U = M x r”, C = M~ l B T , A = BC 


and factors A. Second, one obtains A n by solving 


( 11 ) 


A n = A-\Bu u -rl)/6 2 ( 12 ) 

Finally, u" is obtained from 

u n = u„ - S 2 CX n (13) 

It should be noted that the accuracy loss associated with the factoring of an ill- 
conditioned matrix BA~ 1 B T and the subsequent backsubstitutions can severely influence 
the solution accuracy of not only the Lagrange multipliers but also the generalized coordi- 
nates as seen from (12) and (13). This has motivated many numerical analysts to undertake 
the development of methods for differential-algebraic systems as the recent monograph 46 
and references therein attest to their rich numerical properties. It is generally agreed that 
the present status of differential-algebraic methods yield robust solutions for problems of 
index one , but can suffer from inaccurate solutions of the Lagrange multipliers for higher 
index problems. Observe that many practical multibody dynamics problems are charac- 
terized by index greater than one. Hence, the need to compute accurately the constraint 
forces remains a challenge. For instance, for lock-up mechanisms that are activated when 
truss structures are fully deployed in space often introduce stiff responses with nearly 
singular state of BM 1 B T . It is with these problems for which more robust constraint 
computation algorithms are called for. 

One way to improve the accuracy of constraint force computations is to adopt index 
reduction strategies as discussed in Ref. 46. However, index reduction inevitably intro- 
duces additional system degrees of freedom in the resulting differential-algebraic equations, 
thus destroying the matrix sparsity of (5) in addition to the increased size of the matrix B. 
In what follows we present an alternative approach based on a parabolic regularization of 
the equations for the Lagrange multipliers, which preserves the first row of (5) and enables 
us to solve A from the parabolic differential equations. 


III. Constraint Handling Techniques 

As alluded to in Introduction, techniques for handling the system constraints consti- 
tute a major part of solution procedures for the numerical simulation of multibody dy- 
namics systems. In this section, we will first review the coordinate partitioning technique, 
Baumgarte s technique and the penalty technique. The staggered stabilization procedure 
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which we advocate will then be described in detail. A distinct feature of the staggered 
stabilization procedure is that it can be implemented in a stand-alone module, thus can 
be interfaced not only with the equation solver for rigid-body systems but with that for 
flexible-body systems as well. 


A. Coordinate Partitioning Technique 

In the coordinate partitioning 28 ’ 33 or singular decomposition technique 20 ’ 30 , one se- 
lects a rank sufficient part of B and partitions it as 


B = [Bi B e J, u=[u; q e J (14) 

where the rank of Bi(m x m) is m and the subscripts ( i , e) refer to internal and external 
variables, respectively. First, we express u; in terms of u e as 

u? = B-\r n x -B e u n e ) (15) 

Since we have . T . 

L 7 J { flj } = 0 ( 16 > 

The first row of (10) reduces to 

(Me + T r MiT)q” = r” (17) 


where 


T = B~ l B t , 


Mi 
0 Mel’ 



(18) 


and 


r e 


= r. 


T r r”. + T T MiB\ 


- l r n x 


(19) 


Once one obtains u”, one can obtain u” from (15) and similarly A from (12). Note that 
even though (17) has a smaller dimension than that of (10a), its left-hand side matrix is in 
general full since T given by (18a) is in general full. Hence, unless T is a constant matrix, 
one must refactor the solution matrix in (17) whenever a new T is formed. 


B. Baumgarte’s Technique 

Baumgarte’s technique 22 ’ 29 is based on the observation that the errors committed 
in computing the constraint conditions (3) or (4) can either be critically damped out or 
exponentially decreased as the integration process continues. Mathematically, this can be 
stated for the configuration constraint equation(3) as 

+ 2 a$ + /3$ = 0 ( 20 ) 

or the motion constraint equation(4) as 

$ + 7 $ = 0 ( 21 ) 
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In terms of the general constraint equation augmentation as given by (5b), the pre- 
ceding stabilization is equivalent to modifying c in (5b) accordingly. Hence, the technique 
can be implemented within the standard augmented form of the equations of motion (5) 
However, if BM 1 B T is ill-conditioned, which can happen since B is in general state- 
dependent, the accuracy of generalized constraint force, A, can be considerably degraded. 
This can occur if any two rows of B are physically similar (i.e., when two members form 
a straight line) or numerically close during three-dimensional orientations. 


C. Penalty Technique 

In the two constraint handling techniques outlined so far, the objective was to satisfy 
the constraint condition 


9 = u (22) 

whose differentiated forms were augmented to the equations of motion. In the penalty 
procedure, one adopts 


A = -<?, 


0 


(23) 


as the basic constraint equations instead of the twice-differentiated form adopted in (5). 

It is noted that the penalty formulation tacitly assumes that there will be violations 
of the constraint condition in actual computations as discussed in Lanczos 51 . If one sub- 
stitutes (23) into the governing equations of motion, the resulting equation becomes 


= Q (24) 

A major drawback of the above penalty procedure is that, once an error is committed 
in computing A, there is no compensation scheme by which the drifting of the numerical 
solution can be corrected. This has led to the development of a staggered stabilized 
procedure as described below. 


D. Staggered Stabilization Procedure 

To illustrate this procedure we will consider the case of nonholonomic constraints. 
Instead of substituting the penalty expression directly into the governing equations of 
motion, first we differentiate (23) once to obtain 


A = 


>*!> 


where we assume the penalty parameter, e, to be constant. 
Second, we obtain for ii from (5a) in the form 


u = M - 1 (Q-B r A) 

and substitute it into (25) to yield 


(25) 


(26) 


e\ + BM - 1 B t A = r A = BM^Q + — 

dt 


( 27 ) 
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Notice that the homogeneous part of the above stabilized equation in terms of the 
generalized constraint forces, A, has the following companion eigenvalue problem: 

(7 + BM -1 B r /e)y = 0 ( 28 ) 

where {7*, Jfc = 1 . . . m} are the eigenvalues of the homogeneous operator for the new 
stabilized constraint equations (27). Since 7* also dictates how the errors in the constraint 
forces will diminish with time, the errors committed in the constraint conditions will decay 
with their corresponding different response time constants. This physically oriented stabi- 
lization property of the present technique is in contrast to that of Baumgarte’s technique 
wherein all the error components diminish according to a single time constant. 

Third, this technique enables one to solve for A from the stabilized differential equa- 
tion (27). Specifically, one now has two coupled equations, one set for the generalized 
coordinates u and the other for the generalized constraint forces A, which are recalled here 
from (5a) and (27) for the case of nonholonomic constraints: 


M 

0 


111 !) 


+ 


BM -1 B 


f u 1 _ f Q 1 

J \*J l r A J 


(29) 


Note that the above coupled equations directly provide the desired differential equations 
for a pair of [ u A J . 

For holonomic constraints, one has several stabilization possibilities. The one we have 
chosen is to integrate the governing equations of motion once to obtain 


ii" = (5M _1 (Q" - B r A") + K 


(30) 


which is substituted into Q - 

A=I(B, + f) (3D 

to yield: 

eA” + 6BM -1 B T A n = B(<SM 1 Q n + h?)+ (32) 

It is observed that, even if BM _1 B T is almost singular, this stabilization tech- 
nique as derived in (27) and (32) would not cause numerical difficulty in computing 
A since the solution iteration matrix becomes (e + £BM *B ) for nonholonomic cases 
_|_ <$ 2 BM -1 B r ) for holonomic cases. It is noted that one must choose e in such 
a way to maintain robust solution when BM -1 B T becomes ill-conditioned by choosing 
e ~ c/|(BM - 1 B t ) -1 | • |BM -1 B r | where c is the solution accuracy desired for A. 

Integration of the above equation by the mid-point implicit rule yields the following 
difference equation: 


(el + f BM _1 B T )A n+1/4 = { (r" + * + rj) + eA n (33) 

A n+l/2 = 2A n+1/4 - A" 
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It has been shown that the staggered stabilized procedure for the solution of the 
constraints offers not only a modular software package to treat the constraints but also 
has been found to yield more robust solutions compared to the techniques proposed by 
Baumgarte as reported in Park and Chiou 35 . In particular, even when BM -1 B r be- 
comes nearly singular, the staggered stabilized procedure (33) gives stable and acceptable 
solutions whereas the constraint forces computed by the Baumgarte’s technique diverge. 


IV. Solution Algorithms for Generalized Coordinates 


In addition to the choice of implicit and explicit formulas, the recognition that the 
equations of motion for multibody systems with constraints are not ordinary differential 
equations (ODEs) (see, e.g., Petzold 27 ) has placed a unique requirement in the selection 
of solution algorithms for multibody dynamics problems. From the user’s viewpoint, one 
has the option of either employing one of the available ODE packages (see Enright 32 for 
existing ODE packages) or building a special solution module. It should be noted that, 
since the integration of angular velocity vector does not lead to angular orientations, one 
must solve a set of kinematical equations to obtain the desired angular orientations. 

In this section we describe an explicit-implicit transient analysis algorithm that ex- 
ploits the special kinematical relationships of the generalized rotational coordinates vs. 
the angular velocity, namely, the Euler parameters 34 . The integration of the translational 
coordinates and the angular velocity is accomplished by the central difference formula. It 
should be mentioned that the use of the central difference formula does impose a stepsize 
restriction due to its stability limit ( u > max h < 2) where m mai is the highest angular veloc- 
ity of the system components for rigid-body systems or the highest frequency of the entire 
flexible members for flexible-body systems. The simplicity of its programming effort and 
robustness of its solution results can often become compelling enough to adopt an explicit 
formula, which is the view taken here. 

In conventional structural dynamics analysis, explicit time integration of the equations 
of motion by the central difference formula involves the following two updates per step: 


I u n+1 / 2 = u n_1/2 + /iu n 
l u n+1 = u" + hu n+1/2 


(34) 


Unfortunately, this simplistic procedure is not directly applicable to the rotational part of 
the equations of motion as uj is not directly integrable, except for some special kinematic 
configurations. This motivates us to partition q into the translational velocity vector, d, 
which is directly integrable and the angular velocity vector, u>, which is not, and treat 
them differently, viz.: 

i = {«} (35) 

The equations of motion (5a) can be partitioned according to the above partitioning: 


'M d Olfdl fQrfl 
. 0 J \ ui / { Qu J 


( 36 ) 
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where 


f \ _ / fj-EMd)- Sj(d, e) - BjX 1 
{ Q„ j “ \ tu - D„(w) - S„(d, e) - B.Z A / 


( 37 ) 


in which the subscripts (d,u>) refer to the translational and the rotational motions, re- 
spectively, f is the external force vector, D is the generalized damping force including the 
centrifugal force, S is the internal force vector including member flexibility, q is the angular 
orientation parameters, B d and Bo, are the partition of the combined gradient matrices 
of the constraint conditions (3) or (4) that are symbolically expressed as 


B = Bn + Bh, 


A = A n + 


(38) 


To effect the body-by-body integration for the rotational degrees of freedom, we par- 
tition u; further into 

Ih=|u. 1 ,ih 2 ,...,w'-J T (39) 


where a;^ is a (3x1) angular acceleration vector for the j-th body, 


(40) 


We now present the update algorithm for both translational and rotational coordi- 


nates. 


A. Update of Translational and Angular Velocity 

First, assume that d n+1/2 and q n+1 / 2 are already computed so that we can compute 

d” +1/2 and u> n+1/2 by (36), namely, 


/ d " +1/2 D J +i + s; + ‘ - BjA" + * 

U“ +1/2 J \ d; + * + s; +t - b;a-+" 


(41) 


Second, we update the translational velocity and the angular velocity vectors at the step 
(n+1) by 

Id'"‘ = d" + hd (42) 

l U) 


n+1 -.n -jn+1/2 

d = d +hd 

= u> n + hu; n+1/2 


Third, we update the translational displacement, d, by 


d n+3/2 _ d «+l/2 + ^ 


n+1 


(43) 


However, the updating of the angular orientation requires somewhat involved computa- 
tions. To this end, we will employ the Euler parameters and update them accordingly. 
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B. Update of Euler Parameters and Angular Velocity 

As mentioned in conjunction with a direct use of (34) for integrating the rotational equa- 
tions of motion, it is necessary for one to introduce a set of generalized coordinates whose 
time rate can be related to the angular velocity. To this end, we employ the four-parameter 
Euler representation of the angular velocity for each body as (see, e.g., Wittenburg 10 ): 

q = A(u>)q, q = |.9o 9i 92 93 J T (44) 


q= 2 


0 —u) T 
u —6) 


that is subject to the constraint: 


q T q = 1 


where 


0 

— U>3 

U>2 


0 

-w 1 

— ^2 

UJI 

0 


LO = [o;! U>2 


W 3 J 


T 


and the nodal-designation superscript is omitted for notational simplicity. 

We adopt the mid-point implicit procedure to integrate the Euler parameters: 


' q n+1 = A(o>" +1 )-q n+1 

qti+i = qii+1/2 + Aq^ 1 

■ q **+3/2 _ 2 q "+ 1 _ q **+i/2 
( q *»+ 3 / 2 )r . q n + 3/2 _ j 


(45) 


(46) 


(47) 


It should be noted that the mid-point implicit update is no more costly than any explicit 
as the solution matrix inversion can be explicitly obtained. 

Finally, once q n + 3 / 2 is computed from (47), it is often required to compute the body- 

T* 

fixed basis vector, b = [hi \>2 bsj in terms of the inertial basis vectors, e = 

rp 

[ e x e 2 e 3 J . These two vectors are related by 


where 



b = Re 


(48) 

2(<Zo + 9i) — 1 

2(9i92 + 9093 ) 

2(91 93 - 9092)" 


2(9i92 — 9093 ) 

2(9o + 92 ) — 1 

2(9293 + 9o9i) 

(49) 

. 2(9i93 + 9092 ) 

2(9293 — 9o9i) 

2(9o + 9l) — 1 . 
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C. Update of d,w,d, q at the (n + 2)-step 


So far we have advanced from the step (n+1) to the step (n+3/2). In other words, we have 
advanced only half of the total step. For the next step, viz, the step (n+2) from (n+3/2), 
we employ the following sequence of computations: 


{£}- 


- m-i/ d 2 + '+ s ; +1 -b2v + m 

- M \D ^ +1 +s ; +1 -BjA n+1 J 


/ d” +S/2 = d" +,/2 + hd n+ ' 

l u J n + 3 / 2 = o > n+1 / 2 + /iu ; n+1 

' d"+ 2 = d" +1 + hd n+3/ 2 

q n+3 / 2 = A(w n+ 3 / 2 )q n+3 / 2 

i q n+ 3/2 _ q n+i + Aq"+3/2 

qn-+*2 __ 2q n + 3 / 2 q"+i 

■ (q" +2 ) r q n+2 = 1 


(50) 

(51) 


(52) 


Note that we do not use d”" 1 " 3 ^ 2 and q n+3 / 2 in advancing from the step (n+3/2) to the 
present step (n+2) in computing d n+2 and q" +2 . Instead, we employ d n+1 and q n+1 , hence 
the name two-stage staggered explicit procedure 3 * . The net result is that, even though we 
take a full step (h instead of h/2), we only advance half the step at a time. In other words, 
we evaluate the acceleration and the angular acceleration vectors twice for each full step. 


V. Implementation 

We will now outline the implementation aspects of the the partitioned MBD solution pro- 
cedure. The procedure is implemented into two separate integration modules: generalized- 
coordinate integrator (CINT) and Lagrange multiplier solver (LINT). The generalized- 
coordinate integrator employs a two-stage modified form of the central difference method 
for updating the singular velocity vector and the mid-point implicit rule for updating the 
angular orientations via the Euler parameters. The Lagrange multipliers solver adopts a 
staggered form of the mid-point implicit method. 


A. Generalized-Coordinate Integrator (CINT) 

The module receives /” = B r A n from LINT and advances the solution of the MBD 
equation (1) from time t n to t n+1 . At each integration step, CINT performs the following 
computations. 


Given: 

Compute: 

Advance: 


p" = (d” 1/2 , d n , u> n-1 / 2 , q") and g" = («»,/£ = B T A n )) 
d and w n by (41) 

d n+1/2 = d- 1/2 + Ad" 
d” +1 = d" + Ad" +1/2 
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(53) 



^n+1/2 _ 1/2 ^ ]yjj n 

q" +1 ' 2 = i(I + fA(u> n+1 / 2 )] • q", A = 1 + £(u, 2 + w 2 + W J) 


(54) 


[ q^ 1 = 2q n+1/2 - q”, (q n+1 ) T • q n+1 = 1 

Output: p n+1 = (d” +1 ^ 2 , d n+1 , a> n+1 / 2 , q" +1 ) 

Module Invoke: Call CINT (p n , g n , h , p n+1 ) 

where h is the stepsize and A(u>) is given by 



' 0 

-uq 


— U>3 

1 

UJi 

0 


— W 2 

2 

U) 2 

— W3 

0 



.^3 

L 02 

-<*>i 

0 


(55) 


and q" +1 / 2 is an intermediate vector and (54c) must be solved to obtain q n+1 so as to to 
satisfy the linear dependency constraint, q^q = 1. 


B. Lagrange Multiplier Solver (LINT) 

This module receives (d, d, u>, q) from CINT and performs the following computations. 
Given: ^n+i/2 _ (d n+1 ^ 2 ) d n+1 / 2 , w n+1 / 2 , q n+1 / 2 , A") 

Compute: B n+1 / 2 , BM -1 B T and r" +1/2 by (3) and (4) 

Advance: 

' A" +1/4 = (el + + i(r; + rj +1/2 )) 

' A" + >/ 2 = 2A" +1/ ' 1 - A” (56) 

f"+l/2 _ ^gn+l/2^T . ^n+1/2 


Output: A” +1/2 , f" +1/2 

Module Invoke: Call LINT (£ n+1 / 2 , h, A n+1/2 , f" +1/2 ) 


C. Two-Stage Explicit-Implicit Staggered Procedure 

In order to evaluate u n+1 , u n+1 must be known. Notice from the preceding section that 
only c5 n+1 / 2 is available. Because inaccurate treatments of the gyroscopic damping and the 
centrifugal force terms can lead quickly to computational instability in computing o> n+1 , 
it is not advisable to obtain u> n+1 by extrapolating with o>” +1 / 2 and w n-1 / 2 . To mitigate 
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this difficulty, we advance only to the next half step, at each CINT and LINT call. This 
is illustrated as follows: 


t = t n 

Call CINT (p n , g", h, p n+1 ) 

Call LINT (£" +1 / 2 , h, A n+1/2 , f" +1/2 ) 
t = t n + hj 2 (n <— n + 1/2) 

Call CINT (p" +1/2 , g" +1/2 , h, p n+3 /2) 
Call LINT (^ n+1 , h, A n+1 , f£ +1 ) 
t — t n h 


Note that 


together with 


g n+./2 = („«+!/*, /a ”+>/ 2) 
p n+l/2 = (d”, d” +1/2 , w", q ” +1 / 2 ) 


:»»+l/2 


provides the necessary input data to compute d ' and in the second call of 

CINT in the above calling sequence. In summary, the present procedure requires two 
function evaluations and two A— solutions per each full step, hence the name “two-stage 
explicit-implicit staggered procedure”. 


VI. Numerical Examples 

The two modules, the generalized coordinate integrator (CINT) and the Lagrange multi- 
pliers solver (LINT), have been implemented in Fortran 77. In solving the following three 
example problems, we have incorporated the constraint conditions through the use of La- 
grange multipliers instead of eliminating the constraints. It is therefore necessary to solve 
the governing equations of motion in a way that satisfies the constraint equations. Hence, 
efficient and accurate solutions of these problems will confirm not only the viability of the 
present integration procedure for the solution of the multibody equations of motion with 
or without constraints but also the constraint stabilization procedure in their combined 
totality. 


A. Plane Three-Link Manipulator 

The first problem tested is a simplified version of the seven-link manipulator deployment 
problem 52 . The three links are initially folded and, for modeling simplicity, between the 
two joints is a coil spring which resists a constant deploying force at the tip of the third 
li nk Also, the left-hand end of the first link is fixed through the same coil spring to the 
wall. These three coil springs are to be locked up once the links are deployed straight. The 
deployment sequence of the manipulator is illustrated in Fig. 1. The time-discretized dif- 
ference equations both for Baumgarte’s technique and the staggered stabilization technique 
have been solved at each time increment by a Newton-type iterative procedure to meet 
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a specified accuracy level. Hence, the performance of the two techniques can be assessed 
by the average number of iterations taken per time increment. This is presented in Fig. 
2 for the accuracy of 10 4 . Notice that the staggered stabilization technique requires on 
the average about 4.5 iterations per step, whereas Baumgarte’s technique requires about 
22 iterations per step. 

Note that Baumgarte’s technique fails to converge for time, t sa 1.1 as manifested in Fig. 2 
because the rows in B become numerically dependent upon one another when the links are 
in a straight configuration. This corroborates the theoretical prediction of non- convergence 
whenever the solution matrix, BM -1 B r , for Baumgarte’s technique (see Eqs.(5b), (20) 
and (21)) becomes singular. On the other hand, the staggered stabilization technique still 
converges within 30 iterations, because it overcomes this singularity difficulty, since A still 
exists, as can be seen from Eqs. (27) and (32). 

It should be noted that, in order to avoid such ill-conditioning, one must differen- 
tiate the constraint equations once or twice more and recast the resulting higher-order 
constraint equations in terms of first-order equations with increased number of equations. 
This process is known as an index reduction strategy 46 . Thus, one must restructure the 
augmented equations of motion (5) with the net result of increased solution variables. 
Other techniques involve singular value decompositions, e.g., as advocated by Fiihrer and 
Leimkuhler 53 . On the other hand, the present staggered stabilization technique overcomes 
the ill-conditioning difficulty without restructuring the governing equations of motion. In- 
stead, the constraint equations are enforced in a separate module by the parabolically 
regularized equations for the Lagrange multipliers as derived in (27) and (32). 

Although not reported here, the same relative performance has been observed for different 
accuracy levels, i.e., for the accuracy of 10~ 5 and 10 -6 . 

From this test problem, we conclude that the staggered stabilization technique yields 
both improved accuracy over and greater computational robustness than the Baumgarte 
technique. In addition, the staggered stabilization technique offers software modularity in 
that the solution of the constraint force, A, can be carried out separately from that of the 
generalized displacement, q. The only data each solution module needs to exchange with 
the other is a set of vectors, plus a common module to generate the gradient matrix of the 
constraints, B. However, one should be cautioned not to extrapolate blindly to complex 
problems the results of the present simple examples. Further judicious experiments are 
needed in applying the present staggered stabilization technique to complex production- 

level problems before it can be adopted for general applications in multibody dynamic 
simulations. 


B. Three-Dimensional Double Pendulum 

The second problem with which we have tested the present procedure is a spatially moving 
double pendulum as shown in Fig. 3. The governing equations of motion become those of 
two separate rigid bars, except they are connected by two spherical joints. From Fig. 3 
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we have the the following quantities: 


= d* — -u>* x z* = 0, i — 1, 2. 
2 

(57) 

— 

M = diag{m 1 , J 1 , m 2 , J 2 } 

(58) 

— 

_ fl |z"x 0 0 

B- [l -Az x x -I -|z 2 x 

(59) 

— 


0 
0 
0 

U?2^3( J 2 ~ J 3 ) 
w 3 w l(^3 — Jl) 

U\UJ2{J 1 — J 2 ) 

ii‘ = { d, <i> }' , d' = [r, y, z] T , a;* = [uq, w 2 , ^3] T ( 61 ) 

A = [Ai, A2, A3, A4, A5, \$] T ( 62 ) 

In the preceding equations, is the vectorial distance from the center of the bar to 
the spherical joint constraints, m and J are the three translational and rotatory inertia 
matrices, z is the skew symmetric matrix formed by the three components of z, x implies 
a vector cross multiplication, and the superscript designates the i-th bar. 

The pendulum is originally positioned in a gravity field with initial horizontal angular 
velocities (w? ) = u{ 2) = 1 ). Figure 4 shows the spatial trajectories of the two mass centers 
as projected on the horizontal surface and on the vertical plane. It is noted that the two 
trajectories form a similar pattern. The constraint forces and angular velocities, although 
not reported herein, exhibit patterns that are analogous in their characteristics for the two 
joints and two mass centers, respectively. 

We have performed convergence studies by using different stepsizes h. Numerical evalua- 
tions indicate, as with the rigid-link problem, that when the stepsize samples more than 
20 per period, the present procedure yields both good accuracy and stability. 

C. Open-Loop Torque for Three-Link Manipulator 

The third problem is a three-link manipulator maneuvering under a specified nonholonomic 
jjp velocity constraint. For this problem, both rigid links and flexible links with four 
beam elements per link have been investigated. The flexible beam was modeled with a 
constant-strain Timoshenko beam element that allows large rotations. The three joints are 
modeled as spherical ones and the Lagrange multipliers have been introduced to enforce 
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the joint constraints and well as the nonholonomic constraint at the manipulator tip. The 
trajectories of the manipulator and the tip velocity specification axe shown in Figs. 5 
and 6. The corresponding joint torques for the rigid and flexible links are also shown in 
Figs. 7 and 8, respectively. Note that even though there exists little difference in the two 
trajectories of the rigid and flexible cases, there are significant differences in the open-loop 
joint torques. These will play an important role in the design of controller for vibration 
suppression in the manipulator arms. 

VII. Discussions 

In this chapter, we have presented a computational procedure for direct integration of the 
multibody dynamical (MBD) equations with constraints. 

Because of its step-advancing nature, the procedure is labeled as a two-stage staggered 
explicit-implicit algorithm: explicit for solving the generalized coordinates (CINT) and 
implicit for Lagrange multipliers to incorporate constraints (LINT). Our numerical exper- 
iments indicate that it is essential to enforce the linear dependency constraint condition 
on the Euler parameters at each integration step. 

Numerical experiments reported herein and additional applications conducted so fax in- 
dicate that the present procedure yields robust solutions if the stepsize gives more than 
twenty samples for the period of the apparent highest response frequency of a given multi- 
body system. Hence, the present procedure appears to have accomplished the following: 

• For closed loop multibody systems and/or problems with complex topology wherein it 
is practically inadvisable to eliminate the constraints, the present procedure facilitates 
a straightforward construction of the governing equations of motion with appropri- 
ate constraints. The generalized coordinates and the system open and closed loop 
Lagrange multipliers can then be solved by the present procedure in a partitioned 
manner. 

• For problems that involve lock-up mechanisms or similar discontinuities, the present 
procedure appears to overcome numerical difficulties encountered in using the Baum- 
garte stabilization. This may be an important impetus for applying the present pro- 
cedure for the simulation of deployment dynamics of space structures. 

• The angular velocity is obtained by an adaptation of the central difference algorithm 
in a two-stage form and the update of angular orientations is based on the Euler pa- 
rameters by adopting the mid-point implicit formula. Both of the integrators conserve 
the system energy, which is important when the multibody simulation package is to be 
interfaced with an active control synthesis module. This is because stability margins 
of active control systems are sensitive to the system damping characteristics either 
physical or numerical. 

• The present MBD solution procedure is implemented into two separate modules: the 
generalized coordinates solver (CINT) and the constraint Lagrange multiplier solver 
(LINT). Hence, the task for interfacing of the present MBD solution modules with 
additional capabilities such as active controller, observer and other analysis and design 
software modules becomes relatively straightforward. Such software architecture is 
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in contrast to most of the existing programming practice wherein several analysis 
capabilities axe embedded into a single monolithic program. 

Applications of the present procedure to flexible multibody systems are currently being 
carried out and preliminary results are quite encouraging. We hope to report on the results 
of flexible-body dynamics as well as on large-scale multibody problems in the near future. 
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Fig. 1 Deployment of Three-Link Remote Manipulator 
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Fig 3 Double Pendulum with Spatial Joints 
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Fig. 5 Crane Tip Trajectory of Rigid and Flexible Members 
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Fig. 6 Crane Tip Velocity of Rigid and Flexible Members 



JOINT’TORQUES (lbs in) 




JOINTTORQUES (lbs in) 


1600 


800 


0 


-800 


-1600 



TIME (sec) 


Fig. 8 Cr:ine Joint Torque (Flexible Members) vs. Time 














